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ABSTRACT 

Precision radial velocity (RV) measurements of the Sun-like dwarf 14 Herculis published by 
Naef et. al (2004), Butler et. al (2006) and Wittenmyer et al (2007) reveal a Jovian planet in a 
1760 day orbit and a trend indicating the second distant object. On the grounds of dynamical 
considerations, we test a hypothesis that the trend can be explained by the presence of an 
additional giant planet. We derive dynamical limits to the orbital parameters of the putative 
outer Jovian companion in an orbit within ~ 13 AU. In this case, the mutual interactions 
between the Jovian planets are important for the long-term stability of the system. The best 
self-consistent and stable Newtonian fit to an edge-on configuration of Jovian planets has 
the outer planet in 9 AU orbit with a moderate eccentricity ~ 0.2 and confined to a zone 
spanned by the low-order mean motion resonances 5:1 and 6:1. This solution lies in a shallow 
minimum of (Xv)^^^ persists over a wide range of the system inclination. Other stable 
configurations within la confidence interval of the best fit are possible for the semi-major 
axis of the outer planet in the range of (6,13) AU and the eccentricity in the range of (0,0.3). 
The orbital inclination cannot yet be determined but when it decreases, both planetary masses 
approach ~ lOmj and for / ~ 30^ the hierarchy of the masses is reversed. 
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1 INTRODUCTION 

An extrasolar planet around 14 Herculis discovered by the Geneva 
planet search team was announced in a conference talk (1998). The 
star was also monitored by other planet-hunting teams . The Jovian 
compa nion in ~ 1700 days orbit was next confirm ed bylButler et al.l 
(l2003h . In their new paper, the discovery team (iNaef et al]|2004h 
published 119 observations revealing a linear RV trend of ~ 3.6 m/s 
per year. The single planet Keplerian solution to these measure- 
ments yielded an abnormally large rms of about 14 m/s. Even if the 
drift was accounted for, the rms of the single planet+drift model 
was ~ 11.3 m/s, significantly larger than the mean observational 
uncertainty (Gm) ~ 7.2 m/s quoted in that paper. 

Following the hypothesis that the linea r trend in the RV 
data m ay indicate other, more distant bodies, in lGozdziewski et aP 
we merged these measurement s with accur ate, 35 obser- 
vations (the mean of Gm ~ 3.1 m/s) bv lButler"et aD (2003), span- 
ning the middle part of the combined observational window. We 
reanalyzed the full data set of 154 measureme nts, covering about 
340 days, ~ 2Pb- Using our GAMP method (iGozdziewski et al.l 
1200 3), i.e., the optimization incorporating stability constraints into 
the fitting algorithm, we searched for stable configurations of pu- 
tative two-planet systems which are consistent with the RV data. 
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This search revealed many stable solutions, in particular two dis- 
tinct and equally good best fits with ac ~ 5.8 AU and ac ^ 9 AU. 
They were localized in the proximity of low-order mean motion 
resonances (MMRs), 3:1 MMR (14 Her«) and 6:1 MMR (14 Her^). 
Apparently, the companions would be well separated, but the sys- 
tem would still be active dynamically. The dynamical maps com- 
puted in the neighborhood of the selected fits revealed a complex 
structure of the phase space. We found a clear relation of chaotic 
motions to strongly unstable, self-disrupting configurations on a 
short time- scale of 10^ periods of the outermost planet. 



After the upgrad e of HIRES s pectro ^aph and improvements 
to the data pipeline, iButler et all (l2006h have re- analyzed their 
spectra and published a revised set of 50 observations of 14 Her. 
Their mean accuracy ~ 2.4 m/s is already very good. Yet, towards 
the end of the observational window, si ngle data points have for - 
mal errors as small as 1 m/s. Recently, IWittenmver et al.l (l2007h 
have also published 35 additional and independent precision mea- 
surements made at the McDonald Observatory. They have a mean 
accuracy of 7.5 m/s. The full publicly available data set is now com- 
prised of 203 observations spanning 4463.2 days and still does sup- 
port the presence of a second planetary object. The single-planet 
solution has an rms ~ 13 m/s, and shows a few meter per second 
excess over the me an Gm ^ 6 m /s, added in quadrature to the stel- 
lar jitter ^ 



3.5 m/s (IButler et al. 2006). Using the Keplerian model 
of the RV lWittenmver et al. ( 2007) have found the best-fit solution 
in the vicinity of the 4:1 mean motion resonance (MMR). How- 
ever, having in mind the results of our earlier dynamical analysis 
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jGozdziewski et al.l2006h . we can suspect that the kinematic model 
is not adequate to properly interpret the RV observations. 

In this paper, we re-analyse the updated high-precision RV 
data to study and extend the kinematic model by Wittenmver et al. 
(l2007h . W e also revise and correct t he conclusions from our previ- 
ous paper JGozdziewski et al.ll2006h that were derived on the basis 
of a smaller and less precise RV data set. The goal of our work is 
to derive dynamical characteristics of the putative planetary system 
and to place limits on its orbital parameters. Assuming a putative 
configuration close to the 4:1 MMR, the outermost orbital period 
would be ~ 18-20 yr. Hence, the currently available data would 
cover only a part of this period and a full characterization of the 
orbits would still require many years of observations. One cannot 
expect that in such a case the orbits can be well determined. Nev- 
ertheless, we show in this work, that the modeling of the RV data 
when merged with a dynamical analysis and mapping of the phase 
space of initial conditions with fast-indicators enable us to derive 
meaningful limits to the orbital parameters. Such an extensive dy- 
namical study of this systems has not yet been done. 

In order to take into account the RV variability induced by 
the stellar activity, the internal errors of the RV measurements are 
rescaled by the jitter added in quadrature: = + ' where 
Gm and Gj are the mean of the internal errors and the adopted 
estimate of the st ellar jitter. The 14 Herculis is a quiet star with 
log^HK = -5 07 (iNaef et al.ll2004l) thus it is reas onable to keep the 
safe estimate of Gj = 3.5 m/s (Butler et al. ''2006^. The mass of the 
parent star varies by ~ 10% in different publications. In this work, 
we adopt the canonical mass IMq, although we also did some cal- 
culations for the mass of O.QMq and I.IMq. 



2 MODELING THE RV DATA - AN OVERVIEW 

The standard formulae by 'Sm 3 (11949') are commonly used to 
model the RV signal. Each planet in the system contributes to the 
reflex motion of the star at time t with: 



V,{t) = 7r[cos(co + v(^)) +^cosco] + Vq, 



(1) 



where K is the semi-amplitude, CO is the argument of pericenter, 
v{t) is the true anomaly (involving implicit dependence on the or- 
bital period P and time of periastron passage Tp), e is the eccen- 
tricity, Vq is the velocity offset. We interpret the primary model 
parameters (^,P,^,co, Tp) in terms of the Kepl erian elements and 
minimal masses related to Jacobi coordinates dLee & Pealell2003l : 
iGozdziewski et al. 200 3). 

In our previous works, we tested and optimized different tools 
used to explore the multi-parameter space of x^. In the case when 
has many local extrema, we found that the hybrid optimization 
provides particular ly good results jGozdziewski & Konackill2004l : 
IGozdziewski & Migaszewski 200^^7^ run o f our code starts the 
genetic algorithm (GA. Icharbonneaul Il995h that basically has a 
global nature and requires only the function and rough pa- 
rameter ranges. GA easily permits for a constrained optimization 
within prescribed param eter bounds or for a penalty term to x^ 
(IGozdziewski et al.l 12006) . Because the best fits found with GAs 
are not very accurate in terms of x^, we refine a number of "fittest 
individuals" in the final "population" with a relatively fast local 
method like the simplex of Melder and Nead (Press et al. 19921). 
The use of simplex algorithm is a matter of choice and in princi- 
ple we could use other fast local methods. However, a code using 
non-gradient methods works with the minimal requirements for the 
user- supplied information, i.e., the model function, usually equal to 



X^ [its inverse 1/x^ is the fitness function required by GAs]. The re- 
peated runs of the hybrid code provide an ensemble of the best-fits 
that helps us to detect local minima of x^, even if they are distant in 
the parameter space. We can also o btain reliable approximation to 
the errors of the best-fit parameters (iBevington & Robinson II2003I) 
within the 1g, 2g and 3g confidence levels of x^ at selected 2-dim 
parameter planes, as corresponding to appropriate increments of 
X^. The hybrid approach will be called the algorithm I from here- 
after. 

However, in this paper we deal with a problem of only a partial 
coverage of the longest orbital period by the data. In such a case x^ 
does not have a well defined minimum or the space of acceptable 
solutions is very "flat" so the confidence levels may cover large 
ranges of the parameters. In such a case, to illustrate the shape of 
X^ in the selected 2-dimensional parameter planes, we use a com- 
plementary approach that relies on systematic scanning of the pa- 
rameter space with the fast Levenberg-Marquardt (L-M) algorithm 
(Press et al.lll992 ). Usually, as the parameter plane for represent- 
ing such scans, we choose the semimajor-axis — eccentricity (a, e) 
plane of the outermost planet. We fix (a, e) and then search for the 
best fit, with the initial conditions selected randomly (but within 
reasonably wide parameter bounds). Then the L-M scheme ensures 
a rapid convergence. This approach is c alled the algorithm II from 
hereafter. In fact, we already used it in IGozdziewski et aP (l2003h . 
The method is CPU-time consuming and may be effectively ap- 
plied in low-dimensional problems but in reward it can provide a 
clear and global picture of the parameter space. 

Moreover, the RV measurements can be polluted by many 
sources of error, like complex systematic instrumental effects, short 
time-series of the observations, irregular sampling due to observing 
conditions, and stellar noise. The problem is even more complex 
when we deal with models of resonant or strongly interacting plan- 
etary systems. It is well known that the kinematic model of the 
RV is not adequate to describe the observations in such instances. 
Instea d, the self-consistent A/^-body Newtonian model should be ap- 
pUed dRivera & LissaueJ200ll : lLaughlin & Chambers'2001). (Note 
that the hybrid optimization can be used to explore the param- 
eter space for both models). Still, the best-fit solutions may be 
related (and ofte n do) to unreaUstic quickly disrupting config ura- 
tions (Ferraz-Me llo et al.ll2005l : lLee et al. 2006; Gozdziewski e taP 
l2006l) . In such a case one should explore the dynamical stability 
of the planetary system in a neighborhood of the best-fit configura- 
tion. One can do that with the help of direct numerical integrations 
or by resolving the dynamical character of orbits in terms of the 
maximal Lyapunov exponent or the diffusion of fundamental fre- 
quencies. Hence, the most general way of modeling the RV data 
should involve an elimination of unstable (for instance, strongly 
chaotic) solutions during the fitting procedure. We described such 
an approach (called GAMP) in our previous works. In particular, 
we made attempts to analyze the old data set of 14 Her as well 
as oth er stars hosting multiple planet systems (iGozdziewski et al.! 
l2006h . 



3 KEPLERIAN FITS 



In ord er to compare our fits with the results of I Wittenmver et all 
120071), we carried out the hybrid search for the two-planet edge- 
on configuration. In the experiment, we limited the orbital peri- 
ods to the range of [1000, 14000] days and the eccentricities to the 
range of [0,0.9]. The choi ce of the upper limit of orbital periods 
followed the conclusions of I Wittenmver et al](l2007l) . According to 
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Figure 1. The best fit solutions to the two-planet model of the RV data of 14 Her in the (P, ey and (P,^)-plane. In the search, a distant companion to the known 
Jovian planet in ~ 1760 days orbit is assumed. Upper plots are for the inner planet, middle plots are for the outer one. The best-fit found with the hybrid code in 
the range of moderate orbital periods of the outer planet, with (Xv)^^^ ^ 1.085, is marked with intersecting lines. Its elements, in terms of the parameter tuples 
of Eq.l, (ii:[m/s],P[days],g,co[deg],rp-ro[days]), are (89.460, 1765.743, 0.360, 22.108, 1375.03) for the inner and (226.2263, 11675.333, 0.39 5, 182.006 
6115 .352 ) for the outer planet, re spectively, with the ve locity offsets (Vn2004,^b2006, Vw2007) = (-129.105, -132.345, -108.384) m/s for the data in ( Naef et al.l 
12004) . dButler et al.l t2006) and ( Wittenmver et al."2007), respectively; the epoch 7b = JD2,450,000. Dots are for la solutions with (Xv)^/^ ^ [1.0848, 1.090) 
and an rms ~ 8.17 m/s, red crosses are for the the best fit solutions with (Xv)^^^ ^ 1.0848. The red curve in (Pc,^c)-plane marks an approximate collision 
line of the orbits determined from flb(l + ^b) = «c(l — ^c). with flb^^b fixed at their best-fit values. Bottom panels are for the scans of (Xv)^^^ the selected 
parameter planes. Contours (also shown in the upper panels, for a reference) mark the limits of (Xv)^^^ corresponding to la, 2a and 3a of the best fit solution. 
Circle mark the position of the best A^-body solution (see Sect. 4). See the text for remaining details. 



the res ults of adaptive-optics imaging bv lLuhman & Javawardhana' 
(l2002h and|Patience et aL C200Z) . the parent star has no stellar com- 
panion beyond 9 — 12 AU, so we try to verify the hypothesis about 
a short-period, low-mass companion 14 Her c up to that distance. 
The RV offsets are different for each of the three instruments and 
are included as free parameters in the model together with the tu- 
ples of (^,P,^,co, Tp) for each companion. Note that we calculate 
the RV offset Vn2004 with respect to t he simple mean of all RV mea- 
surements from the set published bv lNaef etal.l(l2004l) . 

The results of the hybrid search (algorithm I) are illustrated in 
upper panels Fig.[T](note that before plotting, the fits are sorted with 
respect to the smaller semi-major axis). The best-fit found with this 
method has (Xv)^^^ — 1085 and an rms ^ 8.17 m/s. Its param- 
eters are given in the caption to Fig. [T] Clearly, there is no iso- 



lated best-fit solution but rather a huge number of acceptable fits. 
The ratio of the orbital periods within la confidence interval of the 
best fit may be as low as ~ 3. The orbital period in the best fit is 
~ 1200 days and PdPh '^1 ■ T his is roughly consistent with the 
work of IWittenmver et all (l2007h who quote the best-fit configura- 
tion close to the 4:1 MMR and with a similar rms 8.17 m/s. Yet 
the fits are not constrained with respect to P^. Many values between 
[5000, 14000] days are equally good in terms of la confidence in- 
terval of the best-fit solution. Moreover, it means that the system 
may be confined to a zone spanned by many low-order mean mo- 
tion resonances of varying width in the (Pc, ^c) -plane. In such cases, 
the mutual interactions between the planets are important for the 
long-term stability. 

In order to illustrate the shape of (Xv)^^^ more detail, we 
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also applied the algorithm II for the systematic scanning of that 
function in selected parameter planes (see lower panels in Fig.[T]). 
This has lead to even better statistics of the best fits and a very clear 
illustration of the (Xv)^^^-surface. The panels reveal that (Xv)^^^ 
of 2-Kepler model has no minimum in the examined ranges of the 
outermost periods, < 14, 000 days. 

Along a wide valley, the best fits with an increasing have 
also rapidly increasing e^. For P^ ~ 14,000 days, they reach the 
collision zone of the orbits. In this strongly chaotic region a num- 
ber of MMRs of the type p : 1 with p >1 overlap. It justifies a- 
posteriori our choice of the search limit P^ ~ 14, 000 days. Note, 
that the two searches are in an excellent accord that is illustrated 
with the (Xv)^^^ smooth contour levels overlay ed on the results of 
the hybrid search. In that case, the scanning of parameter space 
helps to understand the shape of (Xv)^^^ detail. 

The la scatter of solutions around the best-fit (Fig. [B pro- 
vides realistic error estimates. Although the parameters of the inner 
planet are very well fixed, the elements of the outer Jovian planet 
have much larger formal errors. The kinematic fits do not provide 
an upper bound to the semi-major axis of the outer planet. Both 
methods reveal exactly the same parameters of the inner companion 
(see Fig.[T]). Obviously, certain discrepancy between the outcomes 
of algorithm I and algorithm II is related to a "flat" shape of (Xv) 
in the (^c^c) -plane and to a worse efficiency of simplex in detect- 
ing shallow minima compared to the gradient-like L-M algorithm. 
When (Xv)^^^ ^^^^ hdiW^ a well localized minimum, it is very 
difficult to resolve its behaviour with any non-gradient method. 



4 SELF-CONSISTENT NEWTONIAN FITS 

In order to compare the outcome of kinematic modeling with the 
self-consistent N-ho&y one, at first, we transformed the ensemble 
of K eplerian fits to astro-centric osculating elements d Lee & Peald 
12003: Gozdziewski et al. 2003). These fits became the initial con- 
ditions for the L-M algorithm employing the N-bo&y model of the 
RV. In this experiment, the mass of the star was = 1 Mq . Curi- 
ously, we found that the procedure often converged to the two best- 
fits, both having similar (Xv)^^^ ^ 1.083, an rms ~ 8.15 m/s, and 

~ 7.8 AU and ~ 8.4 AU, respectively. To shed some light on 
the orbital stability of these solutions, we computed the dynamical 
maps in the neig hborhood of the best-fits in terms o f the Spectral 
Number flog>SA/^. lMichtchenko & Ferraz^ Melloll200lh . We chose Uc 
and Cc as the map coordinates, keeping other orbital parameters at 
their nominal values. The results are shown in Fig.|2] The left panel 
is for the unstable fit located in the separatrix of the 5: 1 MMR. The 
right panel is for the stable solution in the center of 9:2 MMR. 
These dynamical maps prove that the mutual interactions between 
planets are significant — even small changes of initial parameters 
modify the shapes of low-order MMRs which overlap already for 

~ 0.2 — 0.3. Besides, the chaotic motions are related to strongly 
unstable configurations (see also Fig.[3]|4]). 

As in the case of the Keplerian model, to resolve the topology 
of (Xv)^^^ detail, we performed similar searches for the best-fit 
configurations using the N-hody model of the RV. Again, the re- 
sults of the hybrid algorithm are consistent with the algorithm II. 
However, to conserve space, we only show the (Xv)^^^-niaps ob- 
tained by scanning the (^c^c) -plane. The results are illustrated in 
the left column panels of Fig. [3] 

These scans are computed for three different and fixed values 
of the inclination of the putative orbital plane of the system, / = 
90°, 60° and 30°, respectively. In the maps, we mark a few contour 



levels of (Xv)^^^ that correspond to la, 2a and 3a-levels of the best 
fit solution denoted with a circle in every panel. 

Overall, the topology of (Xv)^^^ ^^^^ change, nevertheless 
qualitatively certain features are common. Contrary to the Keple- 
rian case, we can detect a shallow minimum of (Xv)^^^ close to 
(9 AU,0.2) and a hint of another minimum localized far over the 
collision line. The white thick lines on the panels mark mass lim- 
its of mc. Clearly, for > 9 AU and for moderate eccentricity, rUc 
approaches the brown-dwarf mass limit and substellar mass 80 mj 
for ~ 12 AU. Thus a brown-dwarf could exist in the system beyond 
10-11 AU depending on the inclination and eccentricity — and the 
fits preclude a low-mass planetary object over the 14 mj level. 

Note, that the best-fit for the nominal edge-on configuration 
with / = 90° turns out to be unstable but its parameters are very 
close to a stable solution given in Table 1; its synthetic curve is 
illustrated in Fig. [5] 

To characterize the stability of the best-fit solutions, we com- 
puted their MEGNO fast indicator dCincotta & Sim6ll200Qh . (F), 
(i.e., an approximation of the maximal Lyapunov exponent) up 
to ^ 2 X 10^ Pc using a fast and accurate symplectic algorithm 
jGozdziewski et al.l 12005). The total integration time 3-4 Myr is 
long enough to detect even weak unstable MMRs. 

The best-fits that passed the MEGNO test (i.e. that are reg- 
ular over the integration time) are marked with black dots in the 
right panels of Fig. [3] For a reference, the a-levels of (Xv)^^^ from 
the respective left-panels are marked in these plots as well. The 
red curves are for the collision line of the orbits computed for 
well constrained and fixed parameters of the inner planet. Over- 
all, the stable solutions reveal the structure of the MMRs that we 
have already seen in the dynamical maps computed for fixed or- 
bital parameters (see Fig. |2l). The semi-major axis of the outer 
planet cannot be yet constrained well and may vary between 6 AU 
and (at least) 13 AU. Moreover there is a rather evident although 
irregular border e^ ~ 0.3 of all solutions, up to 3a. Besides, we 
found that in the relevant range of a^, the mean longitude = 
G5c + fAfc ~ (200°, 300°) within the la interval of the best fit, in- 
creasing along quasi-parabolic curve as the function of ^c. Simul- 
taneously, X\y ~ 352.6° ±1.5°. From the dynamical point of view, 
this provides significant limits on the relative mean phases of the 

planets. 

The results for two-planet 14 Her system in lGozdziewski et aP 

(l2006h are basically in accord with the present work although not 
all conclusions are confirmed. The two isolated minima of (Xv)^^^ 
found in that paper can be identified in this work too. In particu- 
lar it concerns the narrow 3:1 MMR island. The second solution 
close to the 6:1 MMR has also been found close to the present 
best fit solution. Nevertheless, the more extensive current search 
reveals many other stable solutions between 3:1 and 9:1 MMRs 
equally good in terms of la. Thus the results of our previous work 
are not wrong in general but rather incomplete — it appears that 
the search performed was not exhaustive enough to find all accept- 
able solutions. Still, some of the unstable solutions can be modi- 
fied in order to obtain stable fits, possibly with slightly increased 
(Xv)^^^- For instance, the stability maps in Fig.[3]seem to rule out 
the 4:1 MMR in the system what may be related to specific orbital 
phases. After an appropriate change of these phases [still keeping 
(Xv)^^^ at acceptable level] we could find also stable solutions close 
to the 4:1 MMR. In order to do this in a self-consistent manner, 
a GAMP-like method would be necessary. However, as we have 
shown, (Xv)^^^ is very flat in the interesting range of orbital pa- 
rameters. Then, the extensive enough GAMP search would require 
a very significant amount of CPU time, so we decided not to do it. 
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Figure 2. The dynamical maps in terms of the Spectral Number (log SN) in the (^c , ^c) -plane for the two-planet Newtonian models of the 14 Her system. Colors 
used in the log SN map classify the orbits — black indicates quasi-periodic regular configurations while yellow strongly chaotic ones. The maps are computed 
for the following osculating elements of the Jovian planets at the epoch of the first observation ^o, in terms of tuples (m [mj],a [AU],e,a) [deg],f7i/ (to) [deg]): 
the left panel is for (4.975, 2.864, 0.359, 22.281, 330.365) of the inner planet, and (5.618, 8.343, 0.174, 188.543, 72.091) of the outer planet, respectively; the 
right panel is for (4.974, 2.863, 0.359, 22.310, 330.314) of the inner planet, and (4.558, 7.863, 0.173, 187.674, 64.417) of the outer planet, respectively. The 
large circles mark the parameters of the best fits. The thin line marks the collision curve for planets b and c, as determined by flb(l + ^b) = flc(l — ^c)- The 
low-order MMRs of the planets b and c are labeled. The integrations are conducted over ~ lO'^Pc- The resolution is 500 x 120 data points. 



Table 1. Astrocentric, osculating Keplerian parameters of the edge-on, 
best-fit, self-consistent configuration of 14 Her system at the epoch of the 
first observation, to = JD 2459464.5956. The mass of the parent star is 
1.0 Mq. See the text and Fig. [3] for the error estimates . The RV offsets 
^N2004» ^62006^ and Vw2007 label the data sets published in jNaef et al.l2004 
IButler et al.i2006 : Wittenm ver et aL 2007) , r espectively; the of fset of yN2004 
is given with respect to the mean of RV in iNaef et al.ll2004h . fA{ denotes 
the mean anomaly at the epoch to. 



parameter 


planet b 




planet c 


m [mj] 


4.975 




7.679 


a [AU] 


2.864 




9.037 


P [days] 


1766 




9886 


e 


0.359 




0.184 


CO [deg] 


22.230 




189.076 


^ (To) [deg] 


330.421 




81.976 


^N2004 [ms-^] 




-52.049 




^B2006 [ms-1] 




-55.290 




Vwiool [ms-^] 




-31.368 








1.0824 




rms [m s ^] 




8.15 





Finally, using the algorithm II, we try to resolve the topology 
of (x?)^^^ with respect to a varying mass of the star and inclina- 
tion of the system (both can be free parameters in the model) but 
assuming that the system remains coplanar. Adopting the uncer- 
tainty of the 14 Her mass ~ 10% dButler et al.ll2006\ we choose 
O.91M0, I.OIMq, I.IMq, and inclinations of orbital plane as 90° 
(the edge-on system), 60° , 45° , 30° . Apparently, the distributions of 
the best-fits are similar for all (M^, /) -pairs. Their quality in terms 
of (Xv)^^^ is comparable. 

The significant differences between the best-fit configurations 
are related to their stability. To address this issue, we calculated 
the dynamical maps for the best fits derived for the nominal mass 
of the parent star (1 Mq) and inclinations 90° (the edge-on sys- 
tem. Table 1), 60°, and 30° (see Fig. Hg]). The log^A^ maps are 
accompanied by the maximal eccentricity maps, max^c, i-e., the 
eccentricity attained by the orbit during the total integration time. 
The max^ indicator makes it possible to resolve a link between 
formally chaotic character of orbits with their physical, short-term 
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Figure 5. Synthetic RV curve of the self-consistent best fit to the two- 
planet model. The osculating elements at the epoch of the first observation 
are given in Table 1. This fit yields (xlY^^ ~ 1.082 and an rms ~ 8.15 m/s. 
Observations from Naef et al. (2004), Butler et al. (2006) and Wittenmyer 
et al. (2007) are marked with red, blue and green circles, accordingly. 

instability. In the regions where chaos is strong, the maximal eccen- 
tricity grows quickly to 1, implying ejections or collisions between 
the planets or with the star. 

The nominal system lies in a stable zone between 5:1 and 
6:1 MMRs (see Fig.|4]). However, for the inclination of ~ 60°, the 
best fit configuration appears locked in the 5:1 MMR. Overall, the 
phase space does not change significantly compared to the edge-on 
system (Fig. ID also Fig. |2]). For the inclination of 30°, the stable 
zone shrinks due to strong interactions caused by large masses of 
the planets ~ lOmj (see caption to Fig. [6]). Moreover, the masses 
are not simply scaled by factor 1 / sin /, as one would expect, assum- 
ing the kinematic model of the RV. Instead, already for / ~ 30°, the 
mass hierarchy is reversed - the inner planet becomes more massive 
than the outer one. It is yet another argument against the kinematic 
model of the RV data. The best-fits derived for small inclinations 
are confined to a strongly chaotic zone. Comparing the dynamical 
map for / ~ 30° with Fig. ID we conclude that is unlikely to find 
stable systems close to the formal best fit with ^ 9 AU unless 
they are trapped in 5:1 MMR. The resonance island can be seen at 
the edge of the bottom panels in Fig. [6] The planetary masses close 
to the minimum of (Xv)^^^ remain smaller from the brown-dwarf 
limit even for relatively small inclinations of the orbital plane. 
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Figure 3. The panels in the left column are for the Newtonian best-fit model with varying fixed inclination of the coplanar system i and illustrated as projections 
onto the (flc,gc) -plane of initial osculating elements at the epoch of the first observation. In each panel, the best fit is marked with a circle. The quality of the 
fits within formal la, 2a and 3a confidence interval of the best fit, in terms of their (Xv)^^^ marked with contours and labeled in subsequent panels. White 
thick lines corresponds to mass levels of a brown-dwarf and sub-stellar companion, respectively. These solutions are derived with algorithm II. The panels in 
the right column are for the stability analysis of the best-fit configurations. See the text for more explanation. 



Still, the parameters of the outer planet are not well con- 
strained and it is clear that many years of new observations are 
required to fix the elements of the putative outer companion. To 
simulate the future RV analysis of the system, we prepared two syn- 
thetic sets of observations, choosing the elements of Table 1 as the 
parameters of the "real" system. We computed the synthetic RV for 
that system and added Gaussian noise to these data with parame- 
ters as from the real observations. Next, we selected the data points 
sampling them with a similar frequency to that in the real measure- 
ments. We analyzed two sets: one extended over -hi 000 days and 
the second one extended over -1-2000 days after the last moment in 
the real observations. Next, for both synthetic data sets (Fig.[7l), we 
scanned the (Xv)^^^ space with the algorithm II in order to resolve 
the best-fit system, moreover, in this experiment, no stability con- 



straints were imposed. The results are shown in two panels in Fig. [8] 
Curiously, after the additional ~ 3 yr observations, the parameters 
of the outer planet cannot be still determined without doubt. After 
~ 6 yr, the space of (Xv) shrinks significantly and we may expect 
that after such a time, the semi-major axis of the outer companion 
can be fixed with la range of ~ 1 AU. 



5 SECULAR DYNAMICS OF THE BEST-FIT SYSTEM 

Although the RV data permit for different dynamical states of the 
system, we can analyse the secular evolution of the best fit coplanar 
and edge-on configuration (Table 1) and the solutions in its neigh- 
borhood. Due to the significant uncertainty of the orbital elements 
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Figure 4. The dynamical maps in the (<3c,^c) -plane in terms of the Spectral Number (log^A^, the top-left panel), maxACH relative to the libration mode 180° 
(the top-right panel) and maxgb,c (the bottom panels) for the two-planet edge-on coplanar Newtonian best-fit model of the 14 Her system (Table 1). Colors used 
in the log SN map classify the orbits — black indicates quasi-periodic regular configurations while yellow strongly chaotic ones. The large crossed circles mark 
the parameters of the fit. The low-order MMRs of the planets b and c are labeled. The integrations are conducted over ~ lO^Pc. The resolution is 500 x 120 
data points. 



of outer planet, such analysis cannot provide exhaustive conclu- 
sions but we can obtain some qualitative view of the secular be- 
haviour of the system. It can be easily extended when more RV 
data will be available. 

As we have demonstrated, the parameters of the inner planet 
as well as the eccentricity and orbital phase of the putative outer 
companion can be relatively well constrained. The most uncon- 
strained parameter in the edge-on, coplanar system is ac. Yet we 
should remember that the observations constrain very roughly the 
inclinations; and we did not include the longitudes of nodes as 
free parameters in the model, keeping the planets in coplanar or- 
bits. Still, we are in relatively a good situation because the care- 
ful analysis of the RV data makes it possible to reduce the di- 
mension of the phase space. The dynamical map in Fig. |4l cover- 
ing the shallow minimum of (Xv)^^^ the neighborhood of the 
formal best fit reveals its proximit y to the 11:2 MMR. Several au- 
thors argue that in similar c ases, (lLibert&Henrardll200l l2005l : 
[Rodriguez & Gallardd l2005l) in the first approximation the near- 
resonance effects have a small or negligible influence on the secular 
dynamics. Hence, neglecting the effects of the MMRs, we can re- 
cover its basic propertie s with the help of the recent qua si-analytical 



this powerful approach can be found in (iMichtchenko & Malhotral 
l2004l : lMichtchenko et al.ll2006h . 



averaging method by Michtchenko & Malhotr a| (l2004l). It general- 
izes t he classic Laplace-Lagrange ( L-L) theorv (iMurrav & DermottI 
I2OOQ) or a more recent work by iLee & Peald ( 1200 3'). It relies 
on a semi-numerical averaging of the perturbing Hamiltonian and 
makes it possible to avoid any series expansions when recovering 
the secular variations of the orbital elements. Then one obtains a 
very accurate secular model which is valid up to large eccentrici- 
ties by far beyond the limits of L-L classical theory. The details of 



secular dynamics, we follow 
j2004h and describe the system in 
First, we obtain the 



To study the 
IMichtchenko & Ma lhotrd 
terms of the canonical Poincare variables, 
time-average of the canonical elements derived from the numerical 
integration of the full equations of motion with the initial condition 
in Table 1. The mean canonical semi-major axes of the system are 
then (ab) 2.868 AU and (ac) 8.891 AU. The eccentricities are 
(^b) ~ 0.356 and (^c) ~ 0.191. The initial state of the system is 
also characterized by the difference of arguments of periastron, 
AG5 = G5b — G5c ~ 195° with la uncertainty ~ 27° that is estimated 
from the statistics of solutions derived for fixed masses and semi- 
major axes (parameters of the secular model), at their respective 
values in Table 1. 

Next, we compute the energy levels of the averaged Hamilto- 
nian of the system (//sec) (Fig- El left panel) in the characteristic 
plane (^bC0sAG5,^c) where AG3 is set either to 0° or 180°. Let us 
note that AG5 always passes through these values during the orbital 
evolution related to two distinct libration modes of AG5 in the copla- 
nar system (IMichtchenko & Malhotral l2004h . The energy level of 
the nominal system is marked with red line. The black thick lines 
mark the exact positions of libration centers (i.e., periodic solu- 
tions) of AG5 around 0° (mode I) and 180° (mode II), respectively. 
A green part of the libration curve (on the right half-plane of the en- 
ergy diagram) is for the nonlinear secular resonance of the system. 
The black dashed lines mark different levels of the total angular 
momentum (expressed by the Angular Momentum Deficit, AMD). 
The nominal system (marked with filled dots) is found in the re- 
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Figure 6. The dynamical maps in the (flc,^c) -plane in terms of the Spectral Number ([ogSN, the left panels), and max^c (the right panels) for the two-planet 
coplanar Newtonian models of the 14 Her system and different inclinations of the orbital plane. Colors used in the \ogSN map classify the orbits — black 
indicates quasi-periodic regular configurations while yellow strongly chaotic ones. The maps are computed for the following elements of the Jovian planets in 
terms of tuples (m [mj],fl [AU],e,a) {to) [deg]). The top row is for i = 60° and the osculating elements (5.747, 2.864, 0.359, 22.271, 330.378) for the 

inner planet, (7.343, 8.675, 0.171, 190.065, 75.540) for the outer planet and the velocity offsets (-43.458, -46.700, -22.780) [m/s]. The fit has {xlY'^ - 1.082 
and rms - 8.15 m/s. The bottom row is for i = 30° and the elements (9.978, 2.868, 0.357, 22.498, 330.135) of the inner planet, (8.581, 7.974, 0.142, 194.823, 
58.154) of the outer planet and the velocity offsets (-29.237, -32.485, -8.585) m/s, [xlf/^ - 1.080 and rms - 8.13 m/s. The large crossed circles mark 
the parameters of the best fits. The low-order MMRs of the planets b and c are labeled. The SN integrations are conducted over ~ lO'^Pc- The resolution is 
500 X 120 data points. Compare with Fig.|4lfor the edge-on system. 
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Figure 7. Synthetic RV data and the RV curve of the corresponding best fit to the two-planet model. The left panel is for the time series extended over 
1000 days after the last real observation, the right panel is for the synthesized data over 2000 days after the last real data point. Sampling frequency and errors 
of the artificial observations are similar to those of the real data. 



gion of the anti-aligned apsides (close to the thin curve representing 
mode II) with in the range ~ (0.1,0.2). 

To illustrate the secular dynamics in more detail, we com- 
puted the phase diagrams of the system in the characteristic plane 
of (^bC0sAG5,^b sinAG5) (the top-left panel in Fig [9]) and in the 
(^cC0sAG5,^c sin AG5) -plane (the bottom-left panel in Fig. [9]). The 
red curves mark the oscillation around libration mode I, blue curves 
are for libration mode II and black curves mark the circulation of 
AG5. The nominal system (marked with filled circles) lies close to 



the center of mode II with a relatively large amplitude of ~ 0.5. 
The librations of AG5 are confined to small semi-amplitudes ~ 30°. 
The results of the secular theory are in an excellent accord with 
the behavior of the full unaveraged system. Its dynamical maps are 
shown in Fig.lH The dynamical map of maxAG5 (i.e., the maximal 
semi- amplitude of AG5 attained during the integration time) reveals 
a similar semi- amplitude of librations ~ 30° around the center of 
180° (and of course, the presence of mode II). This zone persists 
over wide ranges of (ac,^c)- In fact, almost the entire stable region 
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Figure 8. The topology of (Xv)^^^ the (flc,^c) -plane of the best fits to the synthetic data sets shown in Fig. [7] The left panel is for artificial observations 
spanning additional 1000 days after the date of the last real data point, the right panel is for the additional 2000 days period. Contours mark la, 2a and 
3a-levels of (Xy)^^^- Thick lines (white in the left panel, black in the right panel) are for the mass limits of 14 mj and 80 mj, respectively. Circle marks the 
position of the nominal best-fit in Table 1 (the "real" system). 
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Figure 10. The dynamical map in the (^c , ) -plane in terms of the max Acn 
indicator for edge-on, coplanar system and assemble of fits shown in Fig.jS] 
the top-left panel. The large crossed circles mark the parameters of the best 
fit solution in Table 1. The low-order MMRs of the planets b and c are 
labeled. The integrations are conducted over ~ lO^Pc- 



extending at least up to 11 AU is spanned by configurations in- 
volved in mode II librations. The error bars in plots of Fig. [9] are 
derived from formal estimates of la errors of the parameters in Ta- 
ble 1 (we estimate the error of AG5 ~ 27° ; according with the analy- 
sis of stability, we adopted the uncertainty of ec as 0.15). The errors 
are significant but still, in the statistical sense, the libration modes 
would mostly not change within the error ranges. This behaviour is 
confirmed in Fig. [Tol showing AG5 for each individual best-fit found 
(see the top-left panel of Fig. [3]). Simultaneously, physically sta- 
ble systems are close to quasi-periodic configurations what follows 
from the comparison of the max^ and maxAG5 indicators. Both of 
them, as indicators of geometrical or physical behaviour of the sys- 
tem in short-time scale, appear tightly corelated to the log SN as the 
formal measure of the system stability. 



6 CONCLUSIONS 

According to the results from the adaptive-opt ics imaging 
jLuhman & JavawardhanalEo02l : IPatience et al ]|2Q02|), there is no 
stellar-mass object in the 14 Her neighborhood beyond ~ 12 AU. 
This finding supports the hypothesis that the RV trend may be at- 
tributed to a massive planet c. The orbital period of the putative 



companion to the known Jovian planet b cannot yet be very well 
constrained. Nevertheless the available data already reveal a very 
shallow minimum od (Xv)^^^ the (^c^c) -plane of the initial os- 
culating elements. The minimum persists for reasonable combina- 
tions of the parent star mass and the inclination of the system. Its 
dynamical character strongly depends on that parameter influenc- 
ing the planetary masses. Quite surprisingly, for small inclinations 
the mass hierarchy is reversed and the inner planet becomes more 
massive than its distant companion. Also the positions of the nearby 
MMRs are significantly affected. Depending on the inclination, the 
system may be locked in the low-order 9:2, 5:1 or 6:1 MMR. We 
can also conclude that the kinematic model of the RV is already not 
adequate for the characterization of the system and further analysis 
of the RV would be done best within the self-consistent Newtonian 
model. 

The dynamical analysis enables us to derive some limits on 
the orbits elements and the stability of the system. No stable con- 
figurations are found with a period ratio smaller than ~ 3. This 
limit is shifted towards larger values when the inclination grows. 
For the inclination of 30°, the best fit solutions are localized in 
a strongly chaotic and unstable zone. This constitutes a dynami- 
cally derived argument against a small inclination of the system. It 
is likely larger than 30° -40°. Allowing for some speculation, the 
presence of the best-fits in a robustly stable zone supports the hy- 
pothesis of highly inclined two-planet configuration in the 14 Her 
system. Our attempts to fit one more planet to the RV data that 
could explain the scatter of the data points close to the minima of 
the RV curve (Fig. [5]), did not lead to a meaningful improvement of 
the fit. The third planet would have a weak RV signal at the level 
of a few m/s and a short-period eccentric orbit. Finally, the analysis 
of all available RV data of 14 Her enables us to revise and extend 
some of the conclusions from Gozdziewski et al. (2006) — an iso- 
lated and very well defined best-fit solution cannot yet be found to 
properly describe the orbital configuration of 14 Her. 
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Figure 9. The secular evolution of the best-fit configuration (see Table 1). The left panel is for the energy levels of the averaged secular Hamiltonian 
{^)sec (blue, continuous curves) and the AMD levels (black, dashed lines) on the representative plane of initial conditions. The energy decreases and AMD 
increases with increasing ec. The sign of ^bCOsACH is for initial conditions with A03 = 0° (+) or Acn = 180° (-), respectively. The level curves are computed 
for mb/mc = 0.648 and flb/<2c = 0.3226. The red thick curve marks the energy level of the nominal system, and filled circles indicate actual variation of 
the eccentricities (through the angular momentum integral). The black thick lines mark the centers of librations of A03 around 0° (mode I) or 180° (mode II), 
respectively. A green region (on the right, for positive ^b cos A03) is for the nonlinear secular resonance. Panels in the left column are for the secular phase space 
of the planets b and c at secular energy (//)sec of the nominal system in the (gbCOsA03,eb sin Acn) -plane (the right-top panel) and (gcCOsA03,ec sin Acn) -plane 
(the right-bottom panel), respectively. Two libration modes of Acn about the centers AOJ = 0° (mode I, red curves) and Acn = 180° (mode II, blue curves) are 
labeled. Thick black lines are for the transition between the libration and circulation of A03. Initial mean parameters of the nominal systems are marked with 
filled circles. For a reference, we plot formal error bars derived from la errors of the best fit parameters in Table 1. 
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